Tarea_3

# Librerias a utilizar
library(phyloseq)
library(microbiome)

## Loading required package: ggplot2

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>

## 
## Attaching package: 'microbiome'

## The following object is masked from 'package:ggplot2':
## 
##     alpha

## The following object is masked from 'package:base':
## 
##     transform

library(ggplot2)
library(vegan)

## Loading required package: permute

## Loading required package: lattice

## 
## Attaching package: 'vegan'

## The following object is masked from 'package:microbiome':
## 
##     diversity

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Explorando el objeto phyloseq

data("dietswap", package = "microbiome")
filo_sec <- dietswap
filo_sec

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 130 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 130 taxa by 3 taxonomic ranks ]

View(otu_table(filo_sec))

View(sample_data(filo_sec))

View(tax_table(filo_sec)) # no tiene tax table, pero podemos usar la informacion de rank

rank_names(filo_sec)

## [1] "Phylum" "Family" "Genus"

De acuerdo con el codigo anterior ¿Cuántas muestras y taxones contiene el objeto? - Tiene 222 muestras y 130 taxones

¿Qué variables están disponibles en los metadatos de las muestras? - Son 8, las cuales son: numero de muestra, nombre de la muestra, sexo, nacionalidad, grupo, punto de tiempo, punto de tiempo dentro de grupo, indice de peso.

Curvas de rarefacción

Para realizar las curvas de rarefacción usaremos la función rarecurve de vegan.

Para ejecutar la función es necesario convertir nuestra otu table a un data frame o matriz, y posteriormente invertir la posición de las columnas con la de los renglones.

# generamos un objeto a partir de la otu table
otu_table(filo_sec) -> abundancias

# convertimos el objeto a data frame
data_abundancias <- as.data.frame(abundancias)

# transponemos el data frame para que se invierta la posicion de los renglones con las columnas

t(data_abundancias) -> t_data_abundancias

# Para que la curva de cada sample se distinga la una de la otra le pedi a chat gpt un vector de colores

colores_chat <- c("black", "darkblue", "darkcyan", "darkgoldenrod",
                  "darkgray", "darkgreen", "darkgrey", "darkkhaki",
                  "darkmagenta", "darkolivegreen", "darkorange", "darkorchid",
                  "darkred", "darksalmon", "darkseagreen", "darkslateblue",
                  "darkslategray", "darkslategrey", "darkturquoise", "darkviolet",
                  "midnightblue", "navy","brown","deepskyblue", "skyblue")

# Generamos las graficas utilizando el el data frame transpuesto, variando el rango de renglones (samples) que se estan tomando en cuenta.

rarecurve(t_data_abundancias[1:25, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[26:50, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[51:75, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[76:100, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[100:125, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[126:150, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[151:175, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[176:200, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[201:222, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE)

¿Qué indican estas curvas? - Nos dicen si el numero de muestreos por zona fue suficiente para capturar una riqueza de especies cercana a la realidad.

¿Hay muestras que deberían descartarse por bajo conteo? - Si, pero debido a que realice las graficas para conjuntos de 25 muestras no se sabe exactamente que muestra es, de manera que realice el siguiente codigo para averiguarlo.

# Analizando las graficas podemos notar que en el segmento    que va de 51 a 75 hay una muestra que no llega a su         asintota, de manera que trabajare con ese segmento.

t_data_abundancias[51:75,] -> rango 

# Inicializamos la variable x en 50 para que se sume con i y el contador inicie en 51
x <- 50
for (i in 1:25) {
  sum(rango[i,]) -> tamaño_muestral
  x+i -> muestra
  print(paste0("El tamaño muestral es: ",tamaño_muestral,
                " para la muestra ",muestra))
}

## [1] "El tamaño muestral es: 10819 para la muestra 51"
## [1] "El tamaño muestral es: 10979 para la muestra 52"
## [1] "El tamaño muestral es: 8412 para la muestra 53"
## [1] "El tamaño muestral es: 23529 para la muestra 54"
## [1] "El tamaño muestral es: 16757 para la muestra 55"
## [1] "El tamaño muestral es: 1776 para la muestra 56"
## [1] "El tamaño muestral es: 12207 para la muestra 57"
## [1] "El tamaño muestral es: 10503 para la muestra 58"
## [1] "El tamaño muestral es: 6531 para la muestra 59"
## [1] "El tamaño muestral es: 14325 para la muestra 60"
## [1] "El tamaño muestral es: 13984 para la muestra 61"
## [1] "El tamaño muestral es: 8326 para la muestra 62"
## [1] "El tamaño muestral es: 11732 para la muestra 63"
## [1] "El tamaño muestral es: 9614 para la muestra 64"
## [1] "El tamaño muestral es: 8529 para la muestra 65"
## [1] "El tamaño muestral es: 8869 para la muestra 66"
## [1] "El tamaño muestral es: 10920 para la muestra 67"
## [1] "El tamaño muestral es: 8168 para la muestra 68"
## [1] "El tamaño muestral es: 14956 para la muestra 69"
## [1] "El tamaño muestral es: 9328 para la muestra 70"
## [1] "El tamaño muestral es: 10521 para la muestra 71"
## [1] "El tamaño muestral es: 10391 para la muestra 72"
## [1] "El tamaño muestral es: 10011 para la muestra 73"
## [1] "El tamaño muestral es: 12859 para la muestra 74"
## [1] "El tamaño muestral es: 5841 para la muestra 75"

# Nota: Despues de hacer la funcion me di cuenta de que solo debia activar la opcion para ver los nombres dentro del codigo de la misma grafica, pero ya era muy tarde :(

La muestra que se deberia de descartar es la numero 56, ya que cuenta con 1776 individuos y no llega a su asintota. Esto lo comprobamos en la grafica, ya que es la unica linea cuyo tamaño muestral es menor a 5000.

Diversidad alfa (𝛼)

Calcular y graficar los siguientes indices - Riqueza - Shannon - Simpson

# primero recortareos a aquellas taxa que tienen menos de 1 de abundancia
filosec.prune <- prune_taxa(taxa_sums(filo_sec) > 1, filo_sec)

# Riqueza
plot_richness(filosec.prune, x = "nationality", measures = "Observed", color = "sex")

# La grafica nos muestra la riqueza de acuerdo a la         nacionalidad y ditingue los sexos por colores.

# Shannon
plot_richness(filosec.prune, x = "nationality", measures = "shannon", color = "sex")

plot_richness(filosec.prune, measures="Shannon", x="sex", color="bmi_group")

log(108) -> max_sh
max_sh

## [1] 4.682131

# Simpson
plot_richness(filosec.prune, x = "nationality", measures = "simpson", color = "sex")

plot_richness(filosec.prune, measures="simpson", x="sex", color="bmi_group")

¿Qué interpretas de estas gráficas? Que hay unformidad de abundancias entre hombres y mujeres, lo que cambia un poco es la presencia de “lean”, ya que hay más en hombres que en mujeres

¿Hay diferencias notorias entre grupos? Unicamente que hay más “lean” en hombres

Filtrado y transformación

Aplica un filtrado para quedarte solo con los géneros más abundantes de acuerdo a un criterio dado

# Primero crearé un objeto que contenga la riqueza basal y haré comparaciones contra otros objetos
filtro_basal <- prune_taxa(taxa_sums(filo_sec) > 1, filo_sec)
estimate_richness(filtro_basal, measures = "Observed") -> Riqueza_i # Riqueza inicial

filtro_abundancias <- prune_taxa(taxa_sums(filo_sec) > 300, filo_sec)
filtro_abundancias # Este objeto tiene las taxa que tienen una abundancia de al menos 300 en todos los sample

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 98 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 98 taxa by 3 taxonomic ranks ]

estimate_richness(filtro_abundancias, measures = "Observed") -> Riqueza_f # Riqueza final

# Cálculo del cambio de la riqueza de acuerdo al filtrado por abundancias
Riqueza_i-Riqueza_f -> cambio # La riqueza de A debe ser más grande que la de B porque esta considerando más taxa

(cambio/Riqueza_i)*100 -> porcentaje_cambio
min(porcentaje_cambio) # Esto quiere decir que en Riqueza_f tego a los taxones que representan al menos

## [1] 4.901961

# un 4.9% - 100% = 95.6% de la riqueza en los sample
max(porcentaje_cambio) 

## [1] 12.61261

# 12.6% - 100% = 87.4%

# De manera que los generos que representan de un 87.4% a 95.6% de la riqueza se encuntran en el objeto filtro_abundacias
# Dentro del cual tenemos taxa que tienen al menos 300 de abundancia en total 
filtro_abundancias

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 98 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 98 taxa by 3 taxonomic ranks ]

View(otu_table(filtro_abundancias))

###Diversidad beta Realizar una ordención PCoA utilizando distancia Bray-Curtis

# Primero usamos la funcion distance para obtener las distancias con el metodo Bray-curtis
ditancia_braycurtis <- distance(filo_sec, method = "bray")

# Ahora obtenemos las ordenadas
datos_ordenadas_pcoa <- ordinate(filo_sec, method = "PCoA", distance = ditancia_braycurtis)

# primero hacemos el grafico base generando los puntos y separandolos por color de acuerdo a las variables de los metadatos.
# Aquí mismo se generan los ejes de x y y, seleccionado el punto en el que intersectarán.
# Ademas se añadiran los elipses.

pcoa_plot<- plot_ordination(filo_sec, datos_ordenadas_pcoa, 
                             type = "sample", 
                             color = "group") +
  geom_point(size = 0.5) +  
  geom_hline(yintercept = 0, color = "darkgreen", linetype = 5) +
  geom_vline(xintercept = 0, color = "navy", linetype = 5) +
  stat_ellipse(aes(group = group), level = 0.8, color = "darkred", linetype = "dashed") # Los elipses generados son de acuerdo a los elementos dentro de la variable group, es decir : DI, HE y ED.

# Por ultimo se modificaran las etiquetas de los ejes y el titulo de la grafica
pcoa_plot +                         
  ggtitle("PCoA (Bray-Curtis)") +               
  labs(x = paste0("PCoA1 (", round(datos_ordenadas_pcoa$values$Relative_eig[1] * 100, 1), "%)"),
       y = paste0("PCoA2 (", round(datos_ordenadas_pcoa$values$Relative_eig[2] * 100, 1), "%)")) +  # Etiquetas con porcentaje de varianza
  theme_bw() +                                 
  theme(plot.title = element_text(hjust = 0.5)) 

¿Los grupos se separan visiblemente? - no, hay un area muy grande en la que convergen las 3 elipses generadas.

¿Qué podría estar causando esas diferencias? - La ubicación geometrica que tiene cada muestra dentro de las coordenadas principales del grafico. - En este caso las coordenadas principales no son muy explicativas ya que la PCoA 1 explica tan solo un 42% de la varianza, mientras que la PCoA2 explica el 17.2%, eso implica que los grupos DI, HE y ED no son lo suficientemente explicativos como para realizar una agrupación eficiente de las muestras brindadas.

Rank-abundance

# calculo de de la abundancia de cada taxon
abundancias_taxones <- taxa_sums(filo_sec)

# Abundancias ordenadas de mayor a menor
abundancias_ordenadas <- sort(abundancias_taxones, decreasing = TRUE)


# gráfica de barras rank-abundance
par(mgp = c(0, 1, -0.5))
barplot(abundancias_ordenadas, 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "sienna", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7)

¿Qué tan dominada está la comunidad por pocos taxones y cuáles son? Esta muy dominada por pocos taxones, tanto que la mayor abundancia de la muestra se concentra en los primeros 6 taxones de la grafica

¿Qué tipo de distribución parece seguir? logaritmica

Graficas de abundancia por taxón

# Podemos usar únicamente la funcion plot bar
plot_bar(filo_sec, fill = "Phylum")

# o complementarla
plot_bar(filo_sec, fill = "Phylum") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack")

# separada por grupos
plot_bar(filo_sec, fill = "group") + 
  geom_bar(aes(color=group, fill=group), stat="identity", position="stack")

# Agrupe solo a nivel de phylum porque a nivel de genero o familia se ve un desastre.

¿Hay algún phylum que domine? si, los firmicutes.

¿Se observan diferencias entre grupos? Los HE (azules) tienen una mayor abundancia que los ED (verdes), mientras que los DI (rojo claro) son los que tienen menor abundancia.

Exportar resultados

write.csv(abundancias_taxones, file = "../03_Results/suma_abundancias_filosec.csv", row.names = TRUE)

Global patterns

# Preprocesamiento

# Subset para incluir solo muestras de: Soil, Feces, Skin
# Entrega el código y muestra las dimensiones del objeto resultante

data("GlobalPatterns")

GP <- prune_taxa(taxa_sums(GlobalPatterns) > 0, GlobalPatterns)

SFS_subset <- subset_samples(GP, SampleType == c("Soil", "Feces", "Skin"))

## Warning in `==.default`(SampleType, c("Soil", "Feces", "Skin")): longer object
## length is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

SFS_subset

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 18988 taxa and 4 samples ]
## sample_data() Sample Data:       [ 4 samples by 7 sample variables ]
## tax_table()   Taxonomy Table:    [ 18988 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 18988 tips and 18987 internal nodes ]

Diversidad alfa

# Calcular 3 índices de diversidad alfa (Shannon, Simpson, Observed)
indices_diversidad <-plot_richness(GP, x = "Primer", measures = c("Observed", "shannon", "simpson"), color = "SampleType")
indices_diversidad

# Crear boxplots comparativos de los índices entre tipos de muestra
alpha_meas = c("Observed","Shannon", "Simpson")
p <- plot_richness(GP, "human", "SampleType", measures=alpha_meas)

Curvas Rank-abundance

# Crear gráficas de rango-abundancia para cada tipo de muestra
soil_subset <- subset_samples(GP, SampleType == "Soil")
ab_soils <- taxa_sums(soil_subset)
ab_soils_ord  <-sort(ab_soils, decreasing = TRUE)

Feces_subset <- subset_samples(GP, SampleType == "Feces")
ab_feces <- taxa_sums(Feces_subset)
ab_feces_ord <- sort(ab_feces, decreasing = TRUE)

Skin_subset <- subset_samples(GP, SampleType == "Skin")
ab_skin <- taxa_sums(Skin_subset)
ab_skin_ord <- sort(ab_skin, decreasing = TRUE)

Tongue_subset <- subset_samples(GP, SampleType == "Tongue")
ab_tongue <- taxa_sums(Tongue_subset)
ab_tongue_ord <- sort(ab_tongue, decreasing = TRUE)

Freshwater_subset <- subset_samples(GP, SampleType == "Freshwater")
ab_fw <- taxa_sums(Freshwater_subset)
ab_fw_ord <- sort(ab_fw, decreasing = TRUE)

Freshwater_creek_subset <- subset_samples(GP, SampleType == "Freshwater (creek)")
ab_fw_cr <- taxa_sums(Freshwater_creek_subset)
ab_fw_cr_ord <- sort(ab_fw_cr, decreasing = TRUE)

Ocean_subset <- subset_samples(GP, SampleType == "Ocean")
ab_ocean <- taxa_sums(Ocean_subset)
ab_ocean_ord <- sort(ab_ocean, decreasing = TRUE)

Mock_subset <- subset_samples(GP, SampleType == "Mock")
ab_mock <- taxa_sums(Mock_subset)
ab_mock_ord <- sort(ab_mock, decreasing = TRUE)

par(mgp = c(0, 1, -0.5))
barplot(ab_soils_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_feces_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_skin_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_tongue_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_fw_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_fw_cr_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_fw_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_ocean_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

par(mgp = c(0, 1, -0.5))
barplot(ab_mock_ord[1:200], 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "lightseagreen", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7,
        width = 1)

Perfil taxonomico

# Crear gráfico apilado de abundancia a nivel de Phylum
plot_bar(GP, fill = "Phylum") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack")

# Mostrar solo los 5 phyla más abundantes
GP_phylum <- tax_glom(GP, taxrank = "Phylum")
phylum_ab <- taxa_sums(GP_phylum)
cinco <- names(sort(phylum_ab, decreasing = TRUE))[1:5]
seleccion <- prune_taxa(taxa_names(GP_phylum) %in% cinco, GP_phylum)
plot_bar(seleccion, fill = "Phylum") + 
  geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack")

# Agrupar por tipo de muestra
plot_bar(seleccion, fill = "SampleType") + 
  geom_bar(aes(color = SampleType, fill = SampleType), stat = "identity", position = "stack")

# comentar resultados biológicos
# Dado que se realizaron las graficas con los 5 taxones más abundantes pues hay presencia de los 5 taxones en todas las muestras, sin embargo, conforme cambia el tipo de muestra tambien cambia la abundancia de relativa de cada taxón.
# Por ejemplo, en muestras de agua fresca quienes poseen una mayor abundancia son las cianobacterias, mientras que en muestras de tongue quienes tienen mayor abundancia son las proteobacterias

Diversidad beta

                        # Diversidad Beta #

# Calcular distancia Bray-Curtis
d_bc_gp <- distance(GP, method = "bray")

# Ahora obtenemos las ordenadas
datos_pcoa_gp <- ordinate(GP, method = "PCoA", distance = d_bc_gp)

# visualización
pcoa_plot <- plot_ordination(GP, datos_pcoa_gp, 
                             type = "samples",       
                             color = "SampleType") +   
  geom_point(size = 2) +                             
  stat_ellipse(type = "norm", level = 0.95, linetype = 2) +  
  ggtitle("PCoA (Bray-Curtis) de GlobalPatterns") +
  labs(x = paste0("PC1 (", round(datos_pcoa_gp$values$Relative_eig[1]*100, 1), "%)"),
       y = paste0("PC2 (", round(datos_pcoa_gp$values$Relative_eig[2]*100, 1), "%)")) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
pcoa_plot

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_path()`).

#   Incluir stress plot
eig <- datos_pcoa_gp$values$Relative_eig  
scree_plot <- ggplot(data = data.frame(PC = 1:length(eig), Variance = eig),
                     aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  xlab("Eje PCoA") +
  ylab("Proporción de varianza explicada") +
  ggtitle("Scree Plot para PCoA") +
  theme_minimal()
print(scree_plot)

#   Realizar PERMANOVA para diferencias entre grupos
df_gp <- data.frame(sample_data(GlobalPatterns))

# Usamos la función adonis() del paquete vegan:
permanova_result <- adonis(d_bc_gp ~ SampleType, data = df_gp)

## Warning in adonis(d_bc_gp ~ SampleType, data = df_gp): 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").

permanova_result$aov.tab

## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## SampleType  8    7.7744 0.97180  4.4259 0.67562  0.001 ***
## Residuals  17    3.7327 0.21957         0.32438           
## Total      25   11.5070                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Sampletype 8 presenta una diferencia estadisticamente significativa (0.001***) contra el resto de Sampletypes


#   Interpretar resultados en contexto ecológico
# De acuerdon con el PCoA las unicas muestras que se agrupan diferencialmente del resto son las que corresponden a heces (Las 4 que existen son agrpadas en un solo cluster), ya que son las unicas que fueron encerradas en una elipse.
# Esto tambien depende mucho de los ejes, los cuales en este caso explican muy poco de la variación en las muestras (sumados explican menos del 30%), esto se ve reflejado en el scree plot, en el cual vemos que los ejes del 1 al 5 son aquellos que explican entre el 13.9 y el 7.5% de la variación.
#Posteriormente con la aov.tab del permanova observamos que el Sampletye 8 es aquel que muestra una diferencia estadisticamente significativa contra el resto de Sampletypes, lo cual coincide con lo que vemos en el PCoA.

Preguntas de interpretación

¿Qué ambiente muestra mayor diversidad alfa? El suelo.

¿A qué factores podría deberse? A la disponibilidad de nutrientes para los organismos que estan presentes y a condiciones ambientales que favorecen su desarrollo en el mismo.

¿Qué se observa en la ordenación PCoA respecto a similitudes entre muestras? Los unicos grupos que se forman de manera notoria son los de las heces y los de el agua dulce, y esto se debe a la presencia y abundancia de ciertos phylos caracteristicos, tales como las cianobacterias en el aua dulce. Con respecto al resto de grupos que se forman, podemos ver que se agrupan de manera indistinta, ya que hay grupos en los que hay muestras de suelo, oceano y sedimento, sin embargo, la formación de este conjunto no es diferencial.

¿Qué patrón de dominancia taxonómica muestran las curvas de rango-abundancia? Logaritmico o de cola larga

¿Qué phyla son dominantes en cada ambiente y qué implicaciones ecológicas sugieren? agua dulce- arrollo: cianobacterias suelo: proteobacteria mock: Abundancia uniforme de todos lo phyla sedimento: proteobacterias piel: proteobacterias y firmicutes heces: bacteroidetes Tongue: proteobacteria.

Si las muestras son representativas de cada ambiente del que fueron tomadas podria ser un indicador del bacterioma necesario para que la zona de la que fueron tomadas se conserve en el estado en el que estaba cuando se tomaron las muestras.